Handling Geospatial Data in R

Introduction to Geospatial Data in R

This training manual provides a comprehensive guide to handling geospatial data in R, covering both vector and raster data formats, coordinate reference systems, and various mapping techniques.

Learning Objectives

By the end of this tutorial, you will be able to:

• Understand and work with vector data using sf package
• Handle raster data using terra package
• Manage coordinate reference systems (CRS) and projections
• Create static maps with ggplot2 and tmap
• Build interactive maps with leaflet and mapview

Required Packages

First, let’s install and load all necessary packages:

# Install packages if not already installed
# install.packages(c("sf", "terra", "raster", "ggplot2", "tmap", 
#                    "leaflet", "mapview", "dplyr", "viridis", 
#                    "rnaturalearth", "rnaturalearthdata", "lwgeom"))

# Load required libraries
library(sf)              # Simple Features for vector data
library(terra)           # Spatial data analysis (raster)
library(raster)          # Legacy raster package
library(ggplot2)         # Data visualization
library(tmap)            # Thematic mapping
library(leaflet)         # Interactive maps
library(mapview)         # Quick interactive viewing
library(dplyr)           # Data manipulation
library(viridis)         # Color palettes
library(rnaturalearth)   # World map data
library(rnaturalearthdata)

# Set tmap mode to plot (static) initially
tmap_mode("plot")

---

1. Geospatial Vector Data Structure in R

Vector data represents geographic features as points, lines, or polygons with associated attributes.

1.1 Introduction to Simple Features (sf)

The sf package is the modern standard for handling vector geospatial data in R. It integrates seamlessly with the tidyverse and provides a consistent interface.

Key Concepts

• Simple Features: A standardized way to encode spatial vector data
• Geometry Types: POINT, LINESTRING, POLYGON, MULTIPOINT, MULTILINESTRING, MULTIPOLYGON
• sf objects: Data frames with a special geometry column

# Create simple point features manually
points_data <- data.frame(
  name = c("London", "Paris", "Berlin", "Madrid", "Rome"),
  population = c(8982000, 2161000, 3645000, 3223000, 2873000),
  lon = c(-0.1278, 2.3522, 13.4050, -3.7038, 12.4964),
  lat = c(51.5074, 48.8566, 52.5200, 40.4168, 41.9028)
)

# Convert to sf object
cities_sf <- st_as_sf(points_data, 
                      coords = c("lon", "lat"), 
                      crs = 4326)  # WGS84 coordinate system

# Display structure
print(cities_sf)

## Simple feature collection with 5 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -3.7038 ymin: 40.4168 xmax: 13.405 ymax: 52.52
## Geodetic CRS:  WGS 84
##     name population                geometry
## 1 London    8982000 POINT (-0.1278 51.5074)
## 2  Paris    2161000  POINT (2.3522 48.8566)
## 3 Berlin    3645000    POINT (13.405 52.52)
## 4 Madrid    3223000 POINT (-3.7038 40.4168)
## 5   Rome    2873000 POINT (12.4964 41.9028)

# View the geometry column
st_geometry(cities_sf)

## Geometry set for 5 features 
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -3.7038 ymin: 40.4168 xmax: 13.405 ymax: 52.52
## Geodetic CRS:  WGS 84

Exploring sf Object Properties

# Check the geometry type
st_geometry_type(cities_sf)

## [1] POINT POINT POINT POINT POINT
## 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

# Get bounding box
st_bbox(cities_sf)

##    xmin    ymin    xmax    ymax 
## -3.7038 40.4168 13.4050 52.5200

# Get CRS information
st_crs(cities_sf)

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         MEMBER["World Geodetic System 1984 (G2296)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

# Get number of features
nrow(cities_sf)

## [1] 5

# Summary statistics
summary(cities_sf)

##      name             population               geometry
##  Length:5           Min.   :2161000   POINT        :5  
##  Class :character   1st Qu.:2873000   epsg:4326    :0  
##  Mode  :character   Median :3223000   +proj=long...:0  
##                     Mean   :4176800                    
##                     3rd Qu.:3645000                    
##                     Max.   :8982000

1.2 Reading Vector Data from Files

# Get world country boundaries from Natural Earth
world <- ne_countries(scale = "medium", returnclass = "sf")

# Display first few rows
head(world)

## Simple feature collection with 6 features and 168 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.36621 ymin: -22.40205 xmax: 109.4449 ymax: 41.9062
## Geodetic CRS:  WGS 84
##        featurecla scalerank labelrank sovereignt sov_a3 adm0_dif level
## 1 Admin-0 country         1         3   Zimbabwe    ZWE        0     2
## 2 Admin-0 country         1         3     Zambia    ZMB        0     2
## 3 Admin-0 country         1         3      Yemen    YEM        0     2
## 4 Admin-0 country         3         2    Vietnam    VNM        0     2
## 5 Admin-0 country         5         3  Venezuela    VEN        0     2
## 6 Admin-0 country         6         6    Vatican    VAT        0     2
##                type tlc     admin adm0_a3 geou_dif   geounit gu_a3 su_dif
## 1 Sovereign country   1  Zimbabwe     ZWE        0  Zimbabwe   ZWE      0
## 2 Sovereign country   1    Zambia     ZMB        0    Zambia   ZMB      0
## 3 Sovereign country   1     Yemen     YEM        0     Yemen   YEM      0
## 4 Sovereign country   1   Vietnam     VNM        0   Vietnam   VNM      0
## 5 Sovereign country   1 Venezuela     VEN        0 Venezuela   VEN      0
## 6 Sovereign country   1   Vatican     VAT        0   Vatican   VAT      0
##     subunit su_a3 brk_diff      name name_long brk_a3  brk_name brk_group
## 1  Zimbabwe   ZWE        0  Zimbabwe  Zimbabwe    ZWE  Zimbabwe      <NA>
## 2    Zambia   ZMB        0    Zambia    Zambia    ZMB    Zambia      <NA>
## 3     Yemen   YEM        0     Yemen     Yemen    YEM     Yemen      <NA>
## 4   Vietnam   VNM        0   Vietnam   Vietnam    VNM   Vietnam      <NA>
## 5 Venezuela   VEN        0 Venezuela Venezuela    VEN Venezuela      <NA>
## 6   Vatican   VAT        0   Vatican   Vatican    VAT   Vatican      <NA>
##   abbrev postal                        formal_en
## 1  Zimb.     ZW             Republic of Zimbabwe
## 2 Zambia     ZM               Republic of Zambia
## 3   Yem.     YE                Republic of Yemen
## 4  Viet.     VN    Socialist Republic of Vietnam
## 5   Ven.     VE Bolivarian Republic of Venezuela
## 6   Vat.      V        State of the Vatican City
##                            formal_fr              name_ciawf note_adm0 note_brk
## 1                               <NA>                Zimbabwe      <NA>     <NA>
## 2                               <NA>                  Zambia      <NA>     <NA>
## 3                               <NA>                   Yemen      <NA>     <NA>
## 4                               <NA>                 Vietnam      <NA>     <NA>
## 5 República Bolivariana de Venezuela               Venezuela      <NA>     <NA>
## 6                               <NA> Holy See (Vatican City)      <NA>     <NA>
##            name_sort name_alt mapcolor7 mapcolor8 mapcolor9 mapcolor13  pop_est
## 1           Zimbabwe     <NA>         1         5         3          9 14645468
## 2             Zambia     <NA>         5         8         5         13 17861030
## 3        Yemen, Rep.     <NA>         5         3         3         11 29161922
## 4            Vietnam     <NA>         5         6         5          4 96462106
## 5      Venezuela, RB     <NA>         1         3         1          4 28515829
## 6 Vatican (Holy See) Holy See         1         3         4          2      825
##   pop_rank pop_year gdp_md gdp_year                    economy
## 1       14     2019  21440     2019    5. Emerging region: G20
## 2       14     2019  23309     2019  7. Least developed region
## 3       15     2019  22581     2019  7. Least developed region
## 4       16     2019 261921     2019    5. Emerging region: G20
## 5       15     2019 482359     2014    5. Emerging region: G20
## 6        2     2019    -99     2019 2. Developed region: nonG7
##                income_grp fips_10 iso_a2 iso_a2_eh iso_a3 iso_a3_eh iso_n3
## 1           5. Low income      ZI     ZW        ZW    ZWE       ZWE    716
## 2  4. Lower middle income      ZA     ZM        ZM    ZMB       ZMB    894
## 3  4. Lower middle income      YM     YE        YE    YEM       YEM    887
## 4  4. Lower middle income      VM     VN        VN    VNM       VNM    704
## 5  3. Upper middle income      VE     VE        VE    VEN       VEN    862
## 6 2. High income: nonOECD      VT     VA        VA    VAT       VAT    336
##   iso_n3_eh un_a3 wb_a2 wb_a3   woe_id woe_id_eh                   woe_note
## 1       716   716    ZW   ZWE 23425004  23425004 Exact WOE match as country
## 2       894   894    ZM   ZMB 23425003  23425003 Exact WOE match as country
## 3       887   887    RY   YEM 23425002  23425002 Exact WOE match as country
## 4       704   704    VN   VNM 23424984  23424984 Exact WOE match as country
## 5       862   862    VE   VEN 23424982  23424982 Exact WOE match as country
## 6       336   336   -99   -99 23424986  23424986 Exact WOE match as country
##   adm0_iso adm0_diff adm0_tlc adm0_a3_us adm0_a3_fr adm0_a3_ru adm0_a3_es
## 1      ZWE      <NA>      ZWE        ZWE        ZWE        ZWE        ZWE
## 2      ZMB      <NA>      ZMB        ZMB        ZMB        ZMB        ZMB
## 3      YEM      <NA>      YEM        YEM        YEM        YEM        YEM
## 4      VNM      <NA>      VNM        VNM        VNM        VNM        VNM
## 5      VEN      <NA>      VEN        VEN        VEN        VEN        VEN
## 6      VAT      <NA>      VAT        VAT        VAT        VAT        VAT
##   adm0_a3_cn adm0_a3_tw adm0_a3_in adm0_a3_np adm0_a3_pk adm0_a3_de adm0_a3_gb
## 1        ZWE        ZWE        ZWE        ZWE        ZWE        ZWE        ZWE
## 2        ZMB        ZMB        ZMB        ZMB        ZMB        ZMB        ZMB
## 3        YEM        YEM        YEM        YEM        YEM        YEM        YEM
## 4        VNM        VNM        VNM        VNM        VNM        VNM        VNM
## 5        VEN        VEN        VEN        VEN        VEN        VEN        VEN
## 6        VAT        VAT        VAT        VAT        VAT        VAT        VAT
##   adm0_a3_br adm0_a3_il adm0_a3_ps adm0_a3_sa adm0_a3_eg adm0_a3_ma adm0_a3_pt
## 1        ZWE        ZWE        ZWE        ZWE        ZWE        ZWE        ZWE
## 2        ZMB        ZMB        ZMB        ZMB        ZMB        ZMB        ZMB
## 3        YEM        YEM        YEM        YEM        YEM        YEM        YEM
## 4        VNM        VNM        VNM        VNM        VNM        VNM        VNM
## 5        VEN        VEN        VEN        VEN        VEN        VEN        VEN
## 6        VAT        VAT        VAT        VAT        VAT        VAT        VAT
##   adm0_a3_ar adm0_a3_jp adm0_a3_ko adm0_a3_vn adm0_a3_tr adm0_a3_id adm0_a3_pl
## 1        ZWE        ZWE        ZWE        ZWE        ZWE        ZWE        ZWE
## 2        ZMB        ZMB        ZMB        ZMB        ZMB        ZMB        ZMB
## 3        YEM        YEM        YEM        YEM        YEM        YEM        YEM
## 4        VNM        VNM        VNM        VNM        VNM        VNM        VNM
## 5        VEN        VEN        VEN        VEN        VEN        VEN        VEN
## 6        VAT        VAT        VAT        VAT        VAT        VAT        VAT
##   adm0_a3_gr adm0_a3_it adm0_a3_nl adm0_a3_se adm0_a3_bd adm0_a3_ua adm0_a3_un
## 1        ZWE        ZWE        ZWE        ZWE        ZWE        ZWE        -99
## 2        ZMB        ZMB        ZMB        ZMB        ZMB        ZMB        -99
## 3        YEM        YEM        YEM        YEM        YEM        YEM        -99
## 4        VNM        VNM        VNM        VNM        VNM        VNM        -99
## 5        VEN        VEN        VEN        VEN        VEN        VEN        -99
## 6        VAT        VAT        VAT        VAT        VAT        VAT        -99
##   adm0_a3_wb     continent region_un          subregion
## 1        -99        Africa    Africa     Eastern Africa
## 2        -99        Africa    Africa     Eastern Africa
## 3        -99          Asia      Asia       Western Asia
## 4        -99          Asia      Asia South-Eastern Asia
## 5        -99 South America  Americas      South America
## 6        -99        Europe    Europe    Southern Europe
##                    region_wb name_len long_len abbrev_len tiny homepart
## 1         Sub-Saharan Africa        8        8          5  -99        1
## 2         Sub-Saharan Africa        6        6          6  -99        1
## 3 Middle East & North Africa        5        5          4  -99        1
## 4        East Asia & Pacific        7        7          5    2        1
## 5  Latin America & Caribbean        9        9          4  -99        1
## 6      Europe & Central Asia        7        7          4    4        1
##   min_zoom min_label max_label   label_x    label_y      ne_id wikidataid
## 1        0       2.5       8.0  29.92544 -18.911640 1159321441       Q954
## 2        0       3.0       8.0  26.39530 -14.660804 1159321439       Q953
## 3        0       3.0       8.0  45.87438  15.328226 1159321425       Q805
## 4        0       2.0       7.0 105.38729  21.715416 1159321417       Q881
## 5        0       2.5       7.5 -64.59938   7.182476 1159321411       Q717
## 6        0       5.0      10.0  12.45342  41.903323 1159321407       Q237
##     name_ar       name_bn      name_de      name_en             name_es
## 1  زيمبابوي      জিম্বাবুয়ে     Simbabwe     Zimbabwe            Zimbabue
## 2    زامبيا       জাম্বিয়া       Sambia       Zambia              Zambia
## 3     اليمن        ইয়েমেন        Jemen        Yemen               Yemen
## 4    فيتنام      ভিয়েতনাম      Vietnam      Vietnam             Vietnam
## 5   فنزويلا     ভেনেজুয়েলা    Venezuela    Venezuela           Venezuela
## 6 الفاتيكان ভ্যাটিকান সিটি Vatikanstadt Vatican City Ciudad del Vaticano
##    name_fa         name_fr    name_el       name_he   name_hi   name_hu
## 1 زیمبابوه        Zimbabwe Ζιμπάμπουε      זימבבואה   ज़िम्बाब्वे  Zimbabwe
## 2   زامبیا          Zambie     Ζάμπια         זמביה   ज़ाम्बिया    Zambia
## 3      یمن           Yémen     Υεμένη          תימן       यमन     Jemen
## 4   ویتنام        Viêt Nam    Βιετνάμ       וייטנאם   वियतनाम   Vietnám
## 5  ونزوئلا       Venezuela Βενεζουέλα       ונצואלה    वेनेज़ुएला Venezuela
## 6  واتیکان Cité du Vatican   Βατικανό קריית הוותיקן वैटिकन नगर   Vatikán
##     name_id            name_it    name_ja     name_ko      name_nl   name_pl
## 1  Zimbabwe           Zimbabwe ジンバブエ    짐바브웨     Zimbabwe  Zimbabwe
## 2    Zambia             Zambia   ザンビア      잠비아       Zambia    Zambia
## 3     Yaman              Yemen   イエメン        예멘        Jemen     Jemen
## 4   Vietnam            Vietnam   ベトナム      베트남      Vietnam   Wietnam
## 5 Venezuela          Venezuela ベネズエラ  베네수엘라    Venezuela Wenezuela
## 6   Vatikan Città del Vaticano   バチカン 바티칸 시국 Vaticaanstad   Watykan
##     name_pt   name_ru       name_sv   name_tr   name_uk    name_ur
## 1  Zimbábue  Зимбабве      Zimbabwe  Zimbabve  Зімбабве    زمبابوے
## 2    Zâmbia    Замбия        Zambia   Zambiya    Замбія     زیمبیا
## 3     Iémen     Йемен         Jemen     Yemen      Ємен        یمن
## 4  Vietname   Вьетнам       Vietnam   Vietnam   В'єтнам     ویتنام
## 5 Venezuela Венесуэла     Venezuela Venezuela Венесуела  وینیزویلا
## 6  Vaticano   Ватикан Vatikanstaten   Vatikan   Ватикан ویٹیکن سٹی
##         name_vi  name_zh name_zht      fclass_iso tlc_diff      fclass_tlc
## 1      Zimbabwe 津巴布韦   辛巴威 Admin-0 country     <NA> Admin-0 country
## 2        Zambia   赞比亚   尚比亞 Admin-0 country     <NA> Admin-0 country
## 3         Yemen     也门     葉門 Admin-0 country     <NA> Admin-0 country
## 4      Việt Nam     越南     越南 Admin-0 country     <NA> Admin-0 country
## 5     Venezuela 委内瑞拉 委內瑞拉 Admin-0 country     <NA> Admin-0 country
## 6 Thành Vatican   梵蒂冈   梵蒂岡 Admin-0 country     <NA> Admin-0 country
##   fclass_us fclass_fr fclass_ru fclass_es fclass_cn fclass_tw fclass_in
## 1      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 2      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 3      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 4      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 5      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 6      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
##   fclass_np fclass_pk fclass_de fclass_gb fclass_br fclass_il fclass_ps
## 1      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 2      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 3      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 4      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 5      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 6      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
##   fclass_sa fclass_eg fclass_ma fclass_pt fclass_ar fclass_jp fclass_ko
## 1      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 2      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 3      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 4      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 5      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 6      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
##   fclass_vn fclass_tr fclass_id fclass_pl fclass_gr fclass_it fclass_nl
## 1      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 2      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 3      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 4      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 5      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
## 6      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>      <NA>
##   fclass_se fclass_bd fclass_ua                       geometry
## 1      <NA>      <NA>      <NA> MULTIPOLYGON (((31.28789 -2...
## 2      <NA>      <NA>      <NA> MULTIPOLYGON (((30.39609 -1...
## 3      <NA>      <NA>      <NA> MULTIPOLYGON (((53.08564 16...
## 4      <NA>      <NA>      <NA> MULTIPOLYGON (((104.064 10....
## 5      <NA>      <NA>      <NA> MULTIPOLYGON (((-60.82119 9...
## 6      <NA>      <NA>      <NA> MULTIPOLYGON (((12.43916 41...

# Check dimensions
dim(world)

## [1] 242 169

# List column names
names(world)

##   [1] "featurecla" "scalerank"  "labelrank"  "sovereignt" "sov_a3"    
##   [6] "adm0_dif"   "level"      "type"       "tlc"        "admin"     
##  [11] "adm0_a3"    "geou_dif"   "geounit"    "gu_a3"      "su_dif"    
##  [16] "subunit"    "su_a3"      "brk_diff"   "name"       "name_long" 
##  [21] "brk_a3"     "brk_name"   "brk_group"  "abbrev"     "postal"    
##  [26] "formal_en"  "formal_fr"  "name_ciawf" "note_adm0"  "note_brk"  
##  [31] "name_sort"  "name_alt"   "mapcolor7"  "mapcolor8"  "mapcolor9" 
##  [36] "mapcolor13" "pop_est"    "pop_rank"   "pop_year"   "gdp_md"    
##  [41] "gdp_year"   "economy"    "income_grp" "fips_10"    "iso_a2"    
##  [46] "iso_a2_eh"  "iso_a3"     "iso_a3_eh"  "iso_n3"     "iso_n3_eh" 
##  [51] "un_a3"      "wb_a2"      "wb_a3"      "woe_id"     "woe_id_eh" 
##  [56] "woe_note"   "adm0_iso"   "adm0_diff"  "adm0_tlc"   "adm0_a3_us"
##  [61] "adm0_a3_fr" "adm0_a3_ru" "adm0_a3_es" "adm0_a3_cn" "adm0_a3_tw"
##  [66] "adm0_a3_in" "adm0_a3_np" "adm0_a3_pk" "adm0_a3_de" "adm0_a3_gb"
##  [71] "adm0_a3_br" "adm0_a3_il" "adm0_a3_ps" "adm0_a3_sa" "adm0_a3_eg"
##  [76] "adm0_a3_ma" "adm0_a3_pt" "adm0_a3_ar" "adm0_a3_jp" "adm0_a3_ko"
##  [81] "adm0_a3_vn" "adm0_a3_tr" "adm0_a3_id" "adm0_a3_pl" "adm0_a3_gr"
##  [86] "adm0_a3_it" "adm0_a3_nl" "adm0_a3_se" "adm0_a3_bd" "adm0_a3_ua"
##  [91] "adm0_a3_un" "adm0_a3_wb" "continent"  "region_un"  "subregion" 
##  [96] "region_wb"  "name_len"   "long_len"   "abbrev_len" "tiny"      
## [101] "homepart"   "min_zoom"   "min_label"  "max_label"  "label_x"   
## [106] "label_y"    "ne_id"      "wikidataid" "name_ar"    "name_bn"   
## [111] "name_de"    "name_en"    "name_es"    "name_fa"    "name_fr"   
## [116] "name_el"    "name_he"    "name_hi"    "name_hu"    "name_id"   
## [121] "name_it"    "name_ja"    "name_ko"    "name_nl"    "name_pl"   
## [126] "name_pt"    "name_ru"    "name_sv"    "name_tr"    "name_uk"   
## [131] "name_ur"    "name_vi"    "name_zh"    "name_zht"   "fclass_iso"
## [136] "tlc_diff"   "fclass_tlc" "fclass_us"  "fclass_fr"  "fclass_ru" 
## [141] "fclass_es"  "fclass_cn"  "fclass_tw"  "fclass_in"  "fclass_np" 
## [146] "fclass_pk"  "fclass_de"  "fclass_gb"  "fclass_br"  "fclass_il" 
## [151] "fclass_ps"  "fclass_sa"  "fclass_eg"  "fclass_ma"  "fclass_pt" 
## [156] "fclass_ar"  "fclass_jp"  "fclass_ko"  "fclass_vn"  "fclass_tr" 
## [161] "fclass_id"  "fclass_pl"  "fclass_gr"  "fclass_it"  "fclass_nl" 
## [166] "fclass_se"  "fclass_bd"  "fclass_ua"  "geometry"

# Subset to European countries
europe <- world[world$continent == "Europe", ]

# Display basic information
print(paste("Number of European countries:", nrow(europe)))

## [1] "Number of European countries: 50"

print(paste("CRS:", st_crs(europe)$input))

## [1] "CRS: WGS 84"

1.3 Creating Different Geometry Types

Creating Lines

# Create a line connecting cities
city_coords <- st_coordinates(cities_sf)

# Create a LINESTRING
route_matrix <- rbind(
  c(-0.1278, 51.5074),  # London
  c(2.3522, 48.8566),   # Paris
  c(13.4050, 52.5200)   # Berlin
)

route_line <- st_linestring(route_matrix)
route_sf <- st_sf(
  route_name = "London-Paris-Berlin",
  geometry = st_sfc(route_line, crs = 4326)
)

print(route_sf)

## Simple feature collection with 1 feature and 1 field
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -0.1278 ymin: 48.8566 xmax: 13.405 ymax: 52.52
## Geodetic CRS:  WGS 84
##            route_name                       geometry
## 1 London-Paris-Berlin LINESTRING (-0.1278 51.5074...

Creating Polygons

# Create a simple polygon (bounding box around cities)
bbox <- st_bbox(cities_sf)

# Create polygon from bbox
bbox_poly <- st_as_sfc(bbox)
bbox_sf <- st_sf(
  name = "Cities Bounding Box",
  geometry = bbox_poly
)

print(bbox_sf)

## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -3.7038 ymin: 40.4168 xmax: 13.405 ymax: 52.52
## Geodetic CRS:  WGS 84
##                  name                       geometry
## 1 Cities Bounding Box POLYGON ((-3.7038 40.4168, ...

1.4 Spatial Operations with sf

Buffering

# Create 100 km buffer around cities
cities_buffer <- st_buffer(cities_sf, dist = 100000)  # distance in meters

# Plot
plot(st_geometry(europe), col = "lightgray", main = "Cities with 100km Buffers")
plot(st_geometry(cities_buffer), col = "red", alpha = 0.3, add = TRUE)
plot(st_geometry(cities_sf), col = "blue", pch = 19, add = TRUE)
legend("bottomleft", 
       legend = c("Countries", "100km Buffer", "Cities"),
       col = c("lightgray", "red", "blue"),
       pch = c(15, 15, 19))

Spatial Intersection

# Find which countries intersect with city buffers
intersecting <- st_intersection(europe, cities_buffer)

# Count intersections per city
intersection_summary <- intersecting %>%
  group_by(name.1) %>%
  summarise(
    countries_within_100km = n(),
    .groups = 'drop'
  )

print(intersection_summary)

## Simple feature collection with 5 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -4.893563 ymin: 39.50845 xmax: 14.9002 ymax: 53.42907
## Geodetic CRS:  WGS 84
## # A tibble: 5 × 3
##   name.1 countries_within_100km                                         geometry
##   <chr>                   <int>                                    <POLYGON [°]>
## 1 Berlin                      2 ((12.03401 52.18958, 12.02721 52.17413, 12.0415…
## 2 London                      1 ((-0.3889138 52.39264, -0.4011157 52.3926, -0.4…
## 3 Madrid                      1 ((-4.655691 39.88771, -4.655368 39.88772, -4.65…
## 4 Paris                       1 ((2.932405 48.03807, 2.932915 48.04302, 2.96054…
## 5 Rome                        2 ((13.54948 41.46206, 13.55526 41.46137, 13.5552…

Distance Calculations

# Calculate distance matrix between cities (in meters)
dist_matrix <- st_distance(cities_sf)

# Convert to km and create a data frame
dist_km <- units::set_units(dist_matrix, km)
dist_df <- as.data.frame(as.matrix(dist_km))
colnames(dist_df) <- cities_sf$name
rownames(dist_df) <- cities_sf$name

print(round(dist_df, 0))

##           London     Paris    Berlin    Madrid      Rome
## London    0 [km]  344 [km]  932 [km] 1263 [km] 1434 [km]
## Paris   344 [km]    0 [km]  877 [km] 1053 [km] 1105 [km]
## Berlin  932 [km]  877 [km]    0 [km] 1869 [km] 1183 [km]
## Madrid 1263 [km] 1053 [km] 1869 [km]    0 [km] 1364 [km]
## Rome   1434 [km] 1105 [km] 1183 [km] 1364 [km]    0 [km]

Spatial Joins

# Load libraries
library(sf)
library(dplyr)
library(spData)   # spatial example data

# Load country polygons
data("world", package = "spData")

# Load example urban agglomerations (major cities)
data("urban_agglomerations", package = "spData")

# Ensure both layers use the same CRS
urban <- st_transform(urban_agglomerations, st_crs(world))

# Spatial join: city points with country polygons
cities_with_country <- st_join(urban, world, join = st_within)

# Select and rename relevant columns
cities_info <- cities_with_country %>%
  select(
    city = urban_agglomeration,
    country_or_area,
    population_millions,
    continent,        # from world after join
    country = name_long
  ) %>%
  st_drop_geometry()

# Print results
print(cities_info)

## # A tibble: 540 × 5
##    city         country_or_area population_millions continent     country  
##  * <chr>        <chr>                         <dbl> <chr>         <chr>    
##  1 Buenos Aires Argentina                      5.17 South America Argentina
##  2 Buenos Aires Argentina                      5.91 South America Argentina
##  3 Buenos Aires Argentina                      6.76 South America Argentina
##  4 Buenos Aires Argentina                      7.55 South America Argentina
##  5 Buenos Aires Argentina                      8.42 South America Argentina
##  6 Buenos Aires Argentina                      9.14 South America Argentina
##  7 Buenos Aires Argentina                      9.92 South America Argentina
##  8 Buenos Aires Argentina                     10.5  South America Argentina
##  9 Buenos Aires Argentina                     11.1  South America Argentina
## 10 Buenos Aires Argentina                     11.8  South America Argentina
## # ℹ 530 more rows

1.5 Legacy sp Package (Brief Overview)

The sp package was the traditional way to handle spatial data in R. While sf is now preferred, you may encounter sp objects in older code.

# Note: sp is being phased out, but shown for reference
library(sp)

# Convert sf to sp
cities_sp <- as(cities_sf, "Spatial")
class(cities_sp)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

# Convert back to sf
cities_sf_again <- st_as_sf(cities_sp)
class(cities_sf_again)

## [1] "sf"         "data.frame"

---

2. Handling Raster Data in R

Raster data represents continuous spatial phenomena as a grid of cells (pixels), each with a value.

2.1 Introduction to terra Package

The terra package is the modern replacement for the raster package, offering faster performance and improved functionality.

Creating a Simple Raster

# Create a simple raster
r <- rast(nrows = 100, ncols = 100, 
          xmin = -10, xmax = 10, 
          ymin = 40, ymax = 60,
          crs = "EPSG:4326")

# Assign random values
values(r) <- runif(ncell(r), min = 0, max = 100)

# Display raster properties
print(r)

## class       : SpatRaster 
## size        : 100, 100, 1  (nrow, ncol, nlyr)
## resolution  : 0.2, 0.2  (x, y)
## extent      : -10, 10, 40, 60  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :        lyr.1 
## min value   :  0.009598536 
## max value   : 99.997116905

# Plot
plot(r, main = "Random Raster Data", 
     col = viridis(100))

Raster Properties

# Get basic properties
cat("Number of rows:", nrow(r), "\n")

## Number of rows: 100

cat("Number of columns:", ncol(r), "\n")

## Number of columns: 100

cat("Number of cells:", ncell(r), "\n")

## Number of cells: 10000

cat("Resolution:", res(r), "\n")

## Resolution: 0.2 0.2

cat("Extent:", as.vector(ext(r)), "\n")

## Extent: -10 10 40 60

cat("CRS:", crs(r), "\n")

## CRS: GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         MEMBER["World Geodetic System 1984 (G2296)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

# Summary statistics
cat("\nSummary Statistics:\n")

## 
## Summary Statistics:

print(summary(values(r)))

##      lyr.1          
##  Min.   :9.599e-03  
##  1st Qu.:2.450e+01  
##  Median :5.011e+01  
##  Mean   :5.001e+01  
##  3rd Qu.:7.527e+01  
##  Max.   :1.000e+02

2.2 Creating Realistic Raster Data

Elevation-like Surface

# Create an elevation-like surface
x_coords <- seq(-10, 10, length.out = 100)
y_coords <- seq(40, 60, length.out = 100)

# Create a matrix with simulated elevation
elevation_matrix <- outer(x_coords, y_coords, 
                          function(x, y) {
                            100 + 50 * sin(x/2) + 30 * cos(y/3) + 
                            20 * sin(sqrt(x^2 + (y-50)^2)/5)
                          })

# Create raster from matrix
elevation <- rast(elevation_matrix)
ext(elevation) <- c(-10, 10, 40, 60)
crs(elevation) <- "EPSG:4326"
names(elevation) <- "elevation"

# Plot
plot(elevation, 
     main = "Simulated Elevation Surface (m)",
     col = terrain.colors(100))
contour(elevation, add = TRUE, nlevels = 10)

Temperature Surface

# Create a temperature gradient (decreases with latitude)
temp_values <- matrix(nrow = 100, ncol = 100)
for(i in 1:100) {
  for(j in 1:100) {
    lat <- 40 + (60-40) * (i-1)/99
    temp_values[i, j] <- 25 - (lat - 40) * 0.5 + rnorm(1, 0, 1)
  }
}

temperature <- rast(temp_values)
ext(temperature) <- c(-10, 10, 40, 60)
crs(temperature) <- "EPSG:4326"
names(temperature) <- "temperature"

# Plot
plot(temperature, 
     main = "Simulated Temperature in degree Celsius",
     col = rev(heat.colors(100)))

2.3 Raster Operations

Raster Algebra

# Create multiple rasters
r1 <- elevation
r2 <- temperature

# Basic arithmetic
r_sum <- r1 + r2
r_diff <- r1 - r2
r_product <- r1 * r2

# Plot comparison
par(mfrow = c(2, 2))
plot(r1, main = "Elevation", col = terrain.colors(50))
plot(r2, main = "Temperature", col = heat.colors(50))
plot(r_sum, main = "Sum", col = viridis(50))
plot(r_diff, main = "Difference", col = viridis(50))

par(mfrow = c(1, 1))

Raster Reclassification

# Classify elevation into categories
# Low: < 110, Medium: 110-130, High: > 130
rcl_matrix <- matrix(c(
  -Inf, 110, 1,
  110, 130, 2,
  130, Inf, 3
), ncol = 3, byrow = TRUE)

elevation_class <- classify(elevation, rcl_matrix)

# Create labels
levels(elevation_class) <- data.frame(
  value = 1:3,
  category = c("Low", "Medium", "High")
)

# Plot
plot(elevation_class, 
     main = "Elevation Classification",
     col = c("green", "yellow", "brown"))

Focal Operations (Moving Window)

# Calculate mean in 3x3 moving window
elevation_smooth <- focal(elevation, w = 3, fun = mean)

# Compare original and smoothed
par(mfrow = c(1, 2))
plot(elevation, main = "Original Elevation", col = terrain.colors(100))
plot(elevation_smooth, main = "Smoothed Elevation (3x3)", col = terrain.colors(100))

par(mfrow = c(1, 1))

Cropping and Masking

# Define a smaller extent
crop_ext <- ext(-5, 5, 45, 55)

# Crop raster
elevation_crop <- crop(elevation, crop_ext)

# Create a circular mask
center_x <- 0
center_y <- 50
radius <- 5

# Create mask polygon
circle <- st_buffer(st_point(c(center_x, center_y)), dist = radius)
circle_sf <- st_sf(geometry = st_sfc(circle, crs = 4326))

# Mask raster
elevation_mask <- mask(elevation_crop, vect(circle_sf))

# Plot
par(mfrow = c(1, 2))
plot(elevation_crop, main = "Cropped", col = terrain.colors(100))
plot(elevation_mask, main = "Masked (Circular)", col = terrain.colors(100))

par(mfrow = c(1, 1))

2.4 Raster Extraction

# Extract raster values at city locations
cities_elevation <- extract(elevation, vect(cities_sf))
cities_temperature <- extract(temperature, vect(cities_sf))

# Combine with city data
cities_with_raster <- cities_sf %>%
  mutate(
    elevation_m = cities_elevation$elevation,
    temperature_c = cities_temperature$temperature
  ) %>%
  st_drop_geometry()

print(cities_with_raster)

##     name population elevation_m temperature_c
## 1 London    8982000    53.56876      21.18578
## 2  Paris    2161000   141.07056      20.14503
## 3 Berlin    3645000          NA            NA
## 4 Madrid    3223000    39.06862      15.33385
## 5   Rome    2873000          NA            NA

2.5 Raster Stacks and Multi-band Rasters

# Create a multi-layer raster (stack)
environmental_stack <- c(elevation, temperature)
names(environmental_stack) <- c("elevation", "temperature")

print(environmental_stack)

## class       : SpatRaster 
## size        : 100, 100, 2  (nrow, ncol, nlyr)
## resolution  : 0.2, 0.2  (x, y)
## extent      : -10, 10, 40, 60  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## names       : elevation, temperature 
## min values  :  34.38923,    11.76533 
## max values  : 199.86452,    27.60590

# Plot all layers
plot(environmental_stack, 
     col = viridis(100),
     main = c("Elevation (m)", "Temperature in degree Celsius"))

2.6 Raster Package (Legacy)

# The raster package is being phased out
# Converting between terra and raster
library(raster)

# terra to raster
elevation_raster <- raster(elevation)

# raster to terra
elevation_terra <- rast(elevation_raster)

---

3. Projections and Coordinate Reference Systems

Understanding coordinate reference systems (CRS) is crucial for accurate geospatial analysis.

3.1 Introduction to CRS

A Coordinate Reference System defines how spatial data relates to locations on Earth.

Types of CRS

1. Geographic CRS: Uses latitude/longitude (e.g., WGS84, EPSG:4326)
2. Projected CRS: Uses x/y coordinates on a flat surface (e.g., UTM, Web Mercator)

# Check CRS of our data
cat("Cities CRS:\n")

## Cities CRS:

print(st_crs(cities_sf))

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         MEMBER["World Geodetic System 1984 (G2296)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

cat("\nWorld CRS:\n")

## 
## World CRS:

print(st_crs(world))

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

# Check if CRS are the same
cat("\nSame CRS?", st_crs(cities_sf) == st_crs(world), "\n")

## 
## Same CRS? TRUE

3.2 Understanding EPSG Codes

EPSG codes are standardized numerical identifiers for CRS.

Common EPSG Codes

# Create a reference table
epsg_reference <- data.frame(
  EPSG = c(4326, 3857, 27700, 32630, 2154, 3035),
  Name = c("WGS 84", "Web Mercator", "British National Grid", 
           "WGS 84 / UTM zone 30N", "RGF93 / Lambert-93", 
           "ETRS89 / LAEA Europe"),
  Type = c("Geographic", "Projected", "Projected", 
           "Projected", "Projected", "Projected"),
  Use_Case = c("GPS, Global", "Web maps", "UK", 
               "Western Europe", "France", "Europe-wide")
)

print(epsg_reference)

##    EPSG                  Name       Type       Use_Case
## 1  4326                WGS 84 Geographic    GPS, Global
## 2  3857          Web Mercator  Projected       Web maps
## 3 27700 British National Grid  Projected             UK
## 4 32630 WGS 84 / UTM zone 30N  Projected Western Europe
## 5  2154    RGF93 / Lambert-93  Projected         France
## 6  3035  ETRS89 / LAEA Europe  Projected    Europe-wide

Querying CRS Details

# Get detailed information about a CRS
wgs84 <- st_crs(4326)
cat("EPSG 4326 Details:\n")

## EPSG 4326 Details:

cat("Name:", wgs84$Name, "\n")

## Name: WGS 84

cat("Datum:", wgs84$datum, "\n")

## Datum: WGS84

cat("Units:", wgs84$units_gdal, "\n")

## Units: degree

# Check if CRS is geographic or projected
cat("Is geographic?", st_is_longlat(cities_sf), "\n")

## Is geographic? TRUE

3.3 PROJ4 Strings

PROJ4 strings are another way to define CRS, though EPSG codes are preferred.

# Get PROJ4 string
proj4_wgs84 <- st_crs(4326)$proj4string
cat("EPSG:4326 PROJ4 string:\n", proj4_wgs84, "\n\n")

## EPSG:4326 PROJ4 string:
##  +proj=longlat +datum=WGS84 +no_defs

# Create custom projection
custom_crs <- st_crs("+proj=aeqd +lat_0=51.5 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m")
cat("Custom azimuthal equidistant projection:\n")

## Custom azimuthal equidistant projection:

print(custom_crs)

## Coordinate Reference System:
##   User input: +proj=aeqd +lat_0=51.5 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m 
##   wkt:
## PROJCRS["unknown",
##     BASEGEOGCRS["unknown",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]],
##             ID["EPSG",6326]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8901]]],
##     CONVERSION["unknown",
##         METHOD["Azimuthal Equidistant",
##             ID["EPSG",1125]],
##         PARAMETER["Latitude of natural origin",51.5,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["False easting",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]

3.4 Transforming Between CRS

Vector Data Transformation

# Transform cities to different projections
cities_webmercator <- st_transform(cities_sf, crs = 3857)  # Web Mercator
cities_utm <- st_transform(cities_sf, crs = 32630)  # UTM Zone 30N

# Compare coordinates
comparison <- data.frame(
  City = cities_sf$name,
  WGS84_lon = st_coordinates(cities_sf)[,1],
  WGS84_lat = st_coordinates(cities_sf)[,2],
  WebMerc_x = st_coordinates(cities_webmercator)[,1],
  WebMerc_y = st_coordinates(cities_webmercator)[,2],
  UTM_x = st_coordinates(cities_utm)[,1],
  UTM_y = st_coordinates(cities_utm)[,2]
)

print(comparison)

##     City WGS84_lon WGS84_lat  WebMerc_x WebMerc_y     UTM_x   UTM_y
## 1 London   -0.1278   51.5074  -14226.63   6711542  699316.2 5710164
## 2  Paris    2.3522   48.8566  261845.71   6250564  892519.3 5425340
## 3 Berlin   13.4050   52.5200 1492237.77   6894700 1608986.4 5946385
## 4 Madrid   -3.7038   40.4168 -412305.13   4926697  440290.5 4474257
## 5   Rome   12.4964   41.9028 1391092.88   5146430 1786888.1 4756736

Raster Transformation

# Transform elevation raster to UTM
elevation_utm <- project(elevation, "EPSG:32630", method = "bilinear")

# Compare
cat("Original CRS:", crs(elevation, describe = TRUE)$name, "\n")

## Original CRS: WGS 84

cat("Transformed CRS:", crs(elevation_utm, describe = TRUE)$name, "\n")

## Transformed CRS: WGS 84 / UTM zone 30N

# Plot comparison
par(mfrow = c(1, 2))
plot(elevation, main = "WGS84 (Geographic)", col = terrain.colors(100))
plot(elevation_utm, main = "UTM Zone 30N (Projected)", col = terrain.colors(100))

par(mfrow = c(1, 1))

3.5 Area and Distance Calculations

CRS choice affects area and distance calculations significantly.

# Calculate area of European countries in different CRS
europe_wgs84 <- europe
europe_laea <- st_transform(europe, crs = 3035)  # LAEA Europe

# Calculate areas
area_wgs84 <- st_area(europe_wgs84)
area_laea <- st_area(europe_laea)

# Compare for a few countries
comparison_df <- data.frame(
  Country = europe$admin[1:5],
  Area_WGS84_km2 = as.numeric(area_wgs84[1:5]) / 1e6,
  Area_LAEA_km2 = as.numeric(area_laea[1:5]) / 1e6
) %>%
  mutate(Difference_pct = abs(Area_WGS84_km2 - Area_LAEA_km2) / Area_LAEA_km2 * 100)

print(comparison_df)

##          Country Area_WGS84_km2 Area_LAEA_km2 Difference_pct
## 1        Vatican   7.064538e-01  7.075141e-01      0.1498622
## 2         Jersey   1.263326e+02  1.267385e+02      0.3203309
## 3       Guernsey   4.805537e+01  4.821272e+01      0.3263626
## 4    Isle of Man   5.573146e+02  5.597506e+02      0.4351966
## 5 United Kingdom   2.396136e+05  2.406448e+05      0.4285199

cat("\nNote: WGS84 (geographic) gives incorrect areas. Always use projected CRS for area calculations!\n")

## 
## Note: WGS84 (geographic) gives incorrect areas. Always use projected CRS for area calculations!

---

4. Static Maps with ggplot2 and tmap

Creating publication-quality static maps is essential for reports and papers.

4.1 Static Maps with ggplot2 and geom_sf

Basic Map with geom_sf

# Simple map of Europe
ggplot(data = europe) +
  geom_sf(fill = "lightblue", color = "white", size = 0.2) +
  coord_sf(xlim = c(-25, 45), ylim = c(35, 72)) +
  theme_minimal() +
  labs(title = "Map of Europe",
       subtitle = "Using ggplot2 and geom_sf")

Adding Points and Labels

# Map with cities
ggplot() +
  geom_sf(data = europe, fill = "gray90", color = "white", size = 0.3) +
  geom_sf(data = cities_sf, aes(size = population), 
          color = "red", alpha = 0.6) +
  geom_sf_text(data = cities_sf, aes(label = name), 
               nudge_y = 1, size = 3, fontface = "bold") +
  scale_size_continuous(name = "Population", 
                        labels = scales::comma,
                        range = c(2, 10)) +
  coord_sf(xlim = c(-10, 20), ylim = c(40, 60)) +
  theme_minimal() +
  labs(title = "Major European Cities",
       subtitle = "Population size indicated by point size",
       x = "Longitude", y = "Latitude")

Choropleth Map

# Map population density
europe_density <- europe %>%
  mutate(pop_density = pop_est / (as.numeric(st_area(.)) / 1e6))  # per km²

ggplot(data = europe_density) +
  geom_sf(aes(fill = pop_density), color = "white", size = 0.2) +
  scale_fill_viridis_c(name = "Population Density\n(per sq KM)",
                       trans = "log10",
                       labels = scales::comma) +
  coord_sf(xlim = c(-25, 45), ylim = c(35, 72)) +
  theme_minimal() +
  labs(title = "Population Density in Europe",
       subtitle = "Log scale transformation applied")

Multiple Layers

# Complex map with multiple layers
ggplot() +
  geom_sf(data = europe, fill = "gray95", color = "gray70", size = 0.3) +
  geom_sf(data = cities_buffer, fill = "orange", alpha = 0.1, color = "orange") +
  geom_sf(data = route_sf, color = "blue", size = 1.5, linetype = "dashed") +
  geom_sf(data = cities_sf, aes(size = population), 
          color = "red", alpha = 0.8) +
  geom_sf_text(data = cities_sf, aes(label = name), 
               nudge_y = 1.5, size = 3.5, fontface = "bold") +
  scale_size_continuous(name = "Population", 
                        labels = scales::comma,
                        range = c(3, 12)) +
  coord_sf(xlim = c(-12, 18), ylim = c(40, 58)) +
  theme_minimal() +
  theme(legend.position = "bottom") +
  labs(title = "European Cities with Travel Route and Buffers",
       subtitle = "100km buffers shown in orange",
       x = "Longitude", y = "Latitude")

Faceted Maps

# Create categories for population
europe_categories <- europe %>%
  dplyr::filter(continent == "Europe") %>%
  mutate(pop_category = cut(pop_est, 
                            breaks = c(0, 5e6, 20e6, 50e6, Inf),
                            labels = c("< 5M", "5M - 20M", "20M - 50M", "> 50M")))

ggplot(data = europe_categories) +
  geom_sf(aes(fill = pop_category), color = "white", size = 0.2) +
  scale_fill_viridis_d(name = "Population") +
  facet_wrap(~ pop_category) +
  theme_minimal() +
  theme(axis.text = element_blank(),
        axis.ticks = element_blank()) +
  labs(title = "European Countries by Population Category")

4.2 Thematic Mapping with tmap

tmap provides a grammar of graphics specifically designed for maps.

Basic tmap

# Simple map
tm_shape(europe) +
  tm_borders(col = "gray50") +
  tm_fill(col = "lightblue") +
  tm_layout(title = "Europe - Basic tmap")

Choropleth with tmap

# Population map
tm_shape(europe) +
  tm_polygons("pop_est",
              title = "Population",
              palette = "YlOrRd",
              style = "jenks",
              n = 5) +
  tm_layout(title = "European Population",
            legend.outside = TRUE,
            legend.outside.position = "right",
            frame = FALSE)

Advanced tmap with Multiple Elements

# Complex map with multiple layers
tm_shape(europe) +
  tm_polygons("pop_est",
              title = "Population",
              palette = "Blues",
              style = "quantile",
              border.col = "white",
              border.alpha = 0.5) +
  tm_shape(cities_sf) +
  tm_symbols(size = "population",
             col = "red",
             alpha = 0.7,
             title.size = "City Population",
             scale = 1.5) +
  tm_text("name", 
          size = 0.7,
          auto.placement = TRUE,
          fontface = "bold") +
  tm_layout(title = "Europe: Country and City Populations",
            legend.outside = TRUE,
            legend.outside.position = "right",
            bg.color = "lightcyan",
            frame = TRUE) +
  tm_compass(position = c("left", "top")) +
  tm_scale_bar(position = c("left", "bottom"))

Small Multiples with tmap

# Select a few countries for detailed view
selected_countries <- c("France", "Germany", "Spain", "Italy")
countries_subset <- europe %>%
  dplyr::filter(admin %in% selected_countries)

# Faceted map
tm_shape(countries_subset) +
  tm_polygons("pop_est",
              title = "Population",
              palette = "Greens") +
  tm_facets(by = "admin", 
            free.coords = TRUE,
            ncol = 2) +
  tm_layout(legend.show = FALSE,
            panel.labels = selected_countries)

Raster with tmap

# Map elevation raster
tm_shape(elevation) +
  tm_raster(title = "Elevation (m)",
            palette = terrain.colors(100),
            style = "cont") +
  tm_shape(cities_sf) +
  tm_symbols(size = 0.5, col = "red") +
  tm_text("name", size = 0.7, auto.placement = TRUE) +
  tm_layout(title = "Elevation with Cities",
            legend.outside = TRUE,
            bg.color = "lightblue") +
  tm_compass(position = c("right", "top")) +
  tm_scale_bar(position = c("left", "bottom"))

---

5. Interactive Maps with leaflet and mapview

Interactive maps allow users to zoom, pan, and click on features for more information.

5.1 Interactive Maps with leaflet

Basic Leaflet Map

# Simple interactive map
leaflet(data = cities_sf) %>%
  addTiles() %>%  # Add default OpenStreetMap tiles
  addMarkers(popup = ~name)

Customized Markers and Popups

# Create custom popups
popup_content <- paste0(
  "<strong>", cities_sf$name, "</strong><br/>",
  "Population: ", format(cities_sf$population, big.mark = ","), "<br/>",
  "Coordinates: ", round(st_coordinates(cities_sf)[,2], 2), "°N, ",
  round(st_coordinates(cities_sf)[,1], 2), "°E"
)
# Ensure UTF-8
cities_sf$name <- iconv(cities_sf$name, from = "", to = "UTF-8", sub = "")
popup_content <- paste0("<b>", cities_sf$name, "</b><br>Population: ", cities_sf$population)
popup_content <- iconv(popup_content, from = "", to = "UTF-8", sub = "")

# Map with custom popups and circle markers
leaflet(cities_sf) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addCircleMarkers(
    radius = ~sqrt(population) / 100,
    color = "red",
    fillColor = "orange",
    fillOpacity = 0.6,
    popup = popup_content,
    label = ~name
  ) %>%
  addLegend(
    position = "bottomright",
    colors = "orange",
    labels = "European Cities",
    title = "Legend"
  )

Choropleth Map with leaflet

# Create color palette
pal <- colorNumeric(palette = "YlOrRd", domain = europe$pop_est)

# Choropleth map
leaflet(europe) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    fillColor = ~pal(pop_est),
    fillOpacity = 0.7,
    color = "white",
    weight = 1,
    popup = ~paste0("<strong>", admin, "</strong><br/>",
                    "Population: ", format(pop_est, big.mark = ",")),
    highlight = highlightOptions(
      weight = 3,
      color = "red",
      fillOpacity = 0.9,
      bringToFront = TRUE
    )
  ) %>%
  addLegend(
    pal = pal,
    values = ~pop_est,
    title = "Population",
    position = "bottomright",
    labFormat = labelFormat(big.mark = ",")
  )

Multi-layer Interactive Map

# Color palette for raster
raster_pal <- colorNumeric(terrain.colors(100), 
                           values(elevation), 
                           na.color = "transparent")

# Complex map with multiple layers
leaflet() %>%
  addProviderTiles(providers$Esri.WorldImagery, group = "Satellite") %>%
  addProviderTiles(providers$OpenStreetMap, group = "Streets") %>%
  addRasterImage(raster(elevation), 
                 colors = raster_pal, 
                 opacity = 0.6,
                 group = "Elevation") %>%
  addPolygons(data = europe,
              fillColor = "transparent",
              color = "blue",
              weight = 2,
              group = "Countries") %>%
  addCircleMarkers(data = cities_sf,
                   radius = ~sqrt(population) / 100,
                   color = "red",
                   fillColor = "yellow",
                   fillOpacity = 0.8,
                   popup = ~paste0("<strong>", name, "</strong><br/>",
                                   "Pop: ", format(population, big.mark = ",")),
                   group = "Cities") %>%
  addPolylines(data = route_sf,
               color = "purple",
               weight = 3,
               opacity = 0.7,
               group = "Route") %>%
  addLayersControl(
    baseGroups = c("Satellite", "Streets"),
    overlayGroups = c("Elevation", "Countries", "Cities", "Route"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  addLegend(
    pal = raster_pal,
    values = values(elevation),
    title = "Elevation (m)",
    position = "bottomright"
  )

5.2 Quick Interactive Viewing with mapview

mapview provides a simple interface for quick interactive visualization.

Basic mapview

# Quick view of cities
mapview(cities_sf, 
        zcol = "population",
        cex = "population",
        legend = TRUE,
        layer.name = "Population")

Multiple Layers with mapview

# View multiple layers
m1 <- mapview(europe, 
              zcol = "pop_est", 
              layer.name = "Country Pop",
              alpha.regions = 0.5)

m2 <- mapview(cities_sf, 
              zcol = "population",
              col.regions = "red",
              layer.name = "Cities",
              cex = 4)

# Combine maps
m1 + m2

Raster with mapview

# View raster data
mapview(raster(elevation), 
        layer.name = "Elevation",
        col.regions = terrain.colors(100)) +
  mapview(cities_sf, 
          col.regions = "red",
          layer.name = "Cities")

Customized mapview

# More control over appearance
mapview(cities_sf,
        zcol = "population",
        col.regions = c("yellow", "orange", "red"),
        at = c(0, 3000000, 6000000, 9000000),
        cex = 6,
        alpha = 0.8,
        legend = TRUE,
        layer.name = "City Population",
        map.types = c("OpenStreetMap", "Esri.WorldImagery"),
        popup = leafpop::popupTable(
          st_drop_geometry(cities_sf),
          zcol = c("name", "population")
        ))

---

6. Practical Examples and Workflows

6.1 Complete Workflow: Analyzing Accessibility

# Calculate which areas are within 200km of major cities
library(sf)
library(dplyr)
library(ggplot2)
library(viridis)

# --------------------------
# Step 0: Prepare data
# --------------------------
# cities_sf = your 5 major European cities
cities_sf <- st_sf(
  name = c("London", "Paris", "Berlin", "Madrid", "Rome"),
  population = c(8982000, 2161000, 3645000, 3223000, 2873000),
  geometry = st_sfc(
    st_point(c(-0.1278, 51.5074)),
    st_point(c(2.3522, 48.8566)),
    st_point(c(13.405, 52.52)),
    st_point(c(-3.7038, 40.4168)),
    st_point(c(12.4964, 41.9028))
  ),
  crs = 4326
)

# europe = your European country polygons
# europe <- st_read("path_to_europe_shapefile.shp")

# --------------------------
# Step 1: Create 200 km buffers
# --------------------------
cities_200km <- st_transform(cities_sf, 3035) %>%
  st_buffer(200000)

cities_200km_wgs84 <- st_transform(cities_200km, 4326)

# --------------------------
# Step 2: Calculate coverage
# --------------------------
coverage <- st_union(cities_200km_wgs84)

# --------------------------
# Step 3: Find countries with >50% coverage
# --------------------------
europe_3035 <- st_transform(europe, 3035)
cities_200km_3035 <- st_transform(cities_200km_wgs84, 3035)

coverage_analysis <- europe_3035 %>%
  rowwise() %>%
  mutate(
    country_area = as.numeric(st_area(geometry)),
    covered_area = {
      intr <- st_intersection(geometry, st_union(cities_200km_3035))
      intr <- st_collection_extract(intr, "POLYGON")  # only polygons
      if (length(intr) == 0) 0 else as.numeric(st_area(st_union(intr)))
    },
    coverage_pct = (covered_area / country_area) * 100
  ) %>%
  ungroup() %>%
  st_transform(4326) %>%
  arrange(desc(coverage_pct))

# --------------------------
# Step 4: Visualize results
# --------------------------
ggplot() +
  geom_sf(data = coverage_analysis, aes(fill = coverage_pct),
          color = "white", size = 0.2) +
  geom_sf(data = cities_200km_wgs84, fill = "transparent", color = "red",
          linetype = "dashed", alpha = 0.3) +
  geom_sf(data = cities_sf, color = "red", size = 3) +
  scale_fill_viridis_c(name = "Coverage %", option = "plasma") +
  theme_minimal() +
  labs(title = "Accessibility: Areas within 200km of Major Cities",
       subtitle = "Five major European cities analyzed")

# --------------------------
# Step 5: Print top results
# --------------------------
cat("\nTop 10 Countries by Coverage:\n")

## 
## Top 10 Countries by Coverage:

coverage_analysis %>%
  select(admin, coverage_pct) %>%
  st_drop_geometry() %>%
  head(10) %>%
  print()

## # A tibble: 10 × 2
##    admin          coverage_pct
##    <chr>                 <dbl>
##  1 Vatican             100.   
##  2 United Kingdom       29.8  
##  3 Germany              25.3  
##  4 Spain                24.8  
##  5 Italy                20.8  
##  6 France               19.1  
##  7 Poland               10.1  
##  8 Belgium               2.33 
##  9 Czechia               0.678
## 10 Jersey                0

6.2 Raster Analysis Workflow

# Identify high-elevation, low-temperature areas

# Step 1: Define thresholds
high_elevation <- elevation > 130
low_temperature <- temperature < 15

# Step 2: Combine conditions
suitable_areas <- high_elevation & low_temperature

# Step 3: Vectorize suitable areas
suitable_polygons <- as.polygons(suitable_areas) %>%
  st_as_sf() %>%
  dplyr::filter(elevation == 1)  # Keep only TRUE values

# Step 4: Calculate area
suitable_polygons_3035 <- st_transform(suitable_polygons, crs = 3035)
total_area <- sum(as.numeric(st_area(suitable_polygons_3035))) / 1e6  # km²

# Step 5: Visualize
ggplot() +
  geom_sf(data = suitable_polygons, fill = "purple", alpha = 0.5) +
  geom_sf(data = cities_sf, color = "red", size = 3) +
  theme_minimal() +
  labs(title = "Suitable Areas: High Elevation (>130m) & Low Temperature (<15 degree Celsius)",
       subtitle = paste0("Total suitable area: ", round(total_area, 0), " sq KM"))

cat("\nSuitable area statistics:\n")

## 
## Suitable area statistics:

cat("Total area:", round(total_area, 2), "km²\n")

## Total area: 5678.26 km²

cat("Number of polygons:", nrow(suitable_polygons), "\n")

## Number of polygons: 1

---

7. Best Practices and Tips

7.1 Performance Optimization

# 1. Use appropriate CRS for analysis
# - Geographic (EPSG:4326) for global visualization
# - Projected CRS for measurements and analysis

# 2. Simplify geometries when appropriate
simplified <- st_simplify(europe, dTolerance = 1000)

# 3. Use spatial indexing for large datasets
# st_make_valid()  # Fix invalid geometries
# st_is_valid()    # Check validity

# 4. For rasters, use appropriate resolution
# Aggregate to coarser resolution for faster processing
elevation_coarse <- aggregate(elevation, fact = 5, fun = mean)

# 5. Cache intermediate results
# saveRDS(large_sf_object, "cache/processed_data.rds")
# large_sf_object <- readRDS("cache/processed_data.rds")

7.2 Common Pitfalls to Avoid

# 1. NEVER calculate areas in geographic CRS (EPSG:4326)
# WRONG:
area_wrong <- st_area(europe)  # If europe is in EPSG:4326

# RIGHT:
europe_projected <- st_transform(europe, crs = 3035)
area_correct <- st_area(europe_projected)

# 2. Always check CRS compatibility
# st_crs(layer1) == st_crs(layer2)

# 3. Handle invalid geometries
# layer_valid <- st_make_valid(layer_with_issues)

# 4. Be careful with st_intersection on large datasets
# Use st_crop first to reduce computation
# cropped <- st_crop(large_layer, small_extent)
# result <- st_intersection(cropped, other_layer)

# 5. Remember raster vs vector differences
# Rasters: continuous data, fixed resolution
# Vectors: discrete features, variable precision

7.3 Data Export

# Export vector data
st_write(cities_sf, "output/cities.shp")  # Shapefile
st_write(cities_sf, "output/cities.geojson")  # GeoJSON
st_write(cities_sf, "output/cities.gpkg")  # GeoPackage (recommended)

# Export raster data
writeRaster(elevation, "output/elevation.tif")  # GeoTIFF

# Export maps
ggsave("output/map.png", width = 10, height = 8, dpi = 300)
tmap_save(tm_map_object, "output/tmap.png", dpi = 300)

---

Summary and Further Resources

Key Takeaways

1. Vector Data (sf): Modern, tidyverse-compatible package for handling spatial vector data
2. Raster Data (terra): Efficient package for handling gridded spatial data
3. CRS Management: Always use appropriate projections for analysis; EPSG codes are standardized
4. Static Maps: ggplot2 with geom_sf and tmap provide powerful visualization options
5. Interactive Maps: leaflet and mapview enable web-based interactive exploration

Package Summary

Package	Purpose	Key Functions
sf	Vector data	st_read(), st_transform(), st_buffer(), st_intersection()
terra	Raster data	rast(), project(), crop(), mask(), extract()
ggplot2	Static maps	geom_sf(), coord_sf()
tmap	Thematic maps	tm_shape(), tm_polygons(), tm_raster()
leaflet	Interactive maps	leaflet(), addTiles(), addPolygons()
mapview	Quick viewing	mapview()

Additional Resources

• sf documentation: https://r-spatial.github.io/sf/
• terra documentation: https://rspatial.org/terra/
• Geocomputation with R: https://r.geocompx.org/
• Spatial Data Science: https://r-spatial.org/book/
• tmap documentation: https://r-tmap.github.io/tmap/

Practice Exercises

1. Download country boundary data for your region and create a choropleth map
2. Create a 50km buffer around points of interest and calculate coverage statistics
3. Combine multiple raster layers and identify areas meeting specific criteria
4. Build an interactive map with multiple basemaps and layer controls
5. Transform data between different CRS and compare area calculations

---

This material is part of the training program by The National Centre for Research Methods © NCRM authored by Dr Somnath Chaudhuri (University of Southampton). Content is under a CC BY‑style permissive license and can be freely used for educational purposes with proper attribution.